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1 Introduction 

For development of new processes and optimization of existing plants, the 
knowledge of reliable thermophysical data is crucial. Classical methods, next 
to the often expensive and time consuming experimental determination, are 
based on models of Gibbs excess enthalpy (G^-models) or equations of state. 
These phenomenological approaches usually include a large number of ad- 
justable parameters to describe the real behavior quantitatively correct. As in 
most cases these parameters have no physical meaning, their methodic deduc- 
tion is rarely possible and a fit to extensive experimental data is unavoidable. 
Furthermore, extrapolations may be too unreliable for engineering applictions. 

Molecular modeling and simulation offers another approach to thermo- 
physical properties on a much better physical basis. First simulations of this 
type were done by Metropolis et al. in 1953 [1]. By solely defining the interac- 
tions between molecules, it is possible to calculate the relevant thermophysical 
properties with statistical methods. There are two basic techniques which are 
fundamentally different: molecular dynamics (MD) and Monte Carlo (MC). 
The goal of both approaches is to sample the phase space, i.e. all available po- 
sitions and orientations of the molecules under the given boundaries, and an 
estimation of derivatives of the partition function. Macroscopic properties like 
pressure or enthalpy are calculated with statistical thermodynamic methods. 
Detailed descriptions are given, e.g. by Allen and Tildesley [2]. 

In the present project, the molecular interactions are described by effective 
pair potentials. These potentials distinguish dispersive, repulsive, and electro- 
static interactions as well as interactions resulting from hydrogen bonding. 
The dispersive and repulsive part is approximated by Lennard- Jones sites due 
to their computational efficiency. Electrostatics is modeled by partial charges, 
ideal point dipoles and ideal point quadrupoles. For hydrogen bonding good 
results were obtained using eccentric partial charges. 
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Due to a consistent modeling of the pure substances, an application to mix- 
tures is straightforward. Mixed electrostatic potentials are given directly by 
the underlying laws of electrostatics. The parameters of the unlike Lennard- 
Joncs interactions are taken from combining rules, where good results are 
obtained with the Lorentz-Berthelot rule [3, 4]. If needed, one adjustable 
state-independent parameter can be introduced to fit the binary simulation 
results to experimental data. No further parameters were to be used for the 
description of ternary or higher mixtures. 

In the present work, different applications of molecular modeling and sim- 
ulation are presented. Within the project MMSTP, two new molecular models 
for ethylene oxide and ammonia were developed and subsequently used for pre- 
diction of different pure substance properties. Next to this, molecular modeling 
and simulation was applied on the prediction of nucleation processes. Further- 
more, within the present project a comprehensive study on the determination 
of vapor-liquid equilibria of binary and ternary mixtures is performed. This 
is a still ongoing task and thus is not included in the present report. 

Results of this work are consistently published in peer-reviewed interna- 
tional journals. The following publications contribute to the present project: 

• B. Eckl, J. Vrabec, M. Wendland, and H. Hasse: Thermophysical properties 
of dry and humid air by molecular simulation - dew point calculations in 
a new ensemble. In preparation. 

• B. Eckl, J. Vrabec, and H. Hasse: A set of new molecular models based 
on quantum mechanical ab initio calculations. J. Phys. Chem. B, in press 
(2008). 

• M. Horsch, J. Vrabec, and H. Hasse: Molecular dynamics based analysis 
of nucleation and surface energy of droplets in supersaturated vapors of 
methane and ethane. ASME J. Heat Transfer, in press (2008). 

• M. Horsch, J. Vrabec, and H. Hasse: Modification of the classical nucle- 
ation theory based on molecular simulation data for critical nucleus size 
and nucleation rate. Phys. Rev. E, 78: 011603 (2008). 

• M. Horsch, J, Vrabec, M. Bernreuther, S. Grottel, G. Reina, A. Wix, K. 
Schaber, and H. Hasse: Homogeneous nucleation in supersaturated vapors 
of methane, ethane, and carbon dioxide predicted by brute force molecular 
dynamics. J. Chem. Phys., 128: 164510 (2008). 

• B. Eckl, J. Vrabec, and H. Hasse: An optimized molecular model for am- 
monia. Mol. Phys., 106: 1039 (2008). 

• B. Eckl, J. Vrabec, and H. Hasse: On the Application of Force Fields for 
Predicting a Wide Variety of Properties: Ethylene Oxide as an Example. 
Fluid Phase Equilib., in press (2008). 

The present status report briefly shows results of the pure substances 
ethylene oxide and ammonia, demonstrating the outstanding predictive power 
of the applied methods. Secondly, the application to nucleation processes is 
presented. Next to a pure description, the authors were able to further study 
the underlying effects on the molecular level and to deduce a modification to 
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classical nucleation theory. Futher details are given in the references mentioned 
above. 

2 Prediction of a wide variety of properties for ethylene 
oxide 

The predictive and extrapolative power of molecular modeling and simula- 
tion is of particular interest for industrial applications, where a broad variety 
of properties is needed but often not available. To discuss this issue, the In- 
dustrial Simulation Collective [5] has organized the Fourth Industrial Fluid 
Properties Simulation Challenge as an international contest. The task in 2007 
was to calculate for ethylene oxide on the basis of a single molecular model a 
total of 17 different properties from three categories. It should be noted that 
our contribution was awarded with the first price [5] . 

Ethylene oxide (C2H4O) is a widely used intermediate in the chemical in- 
dustry. In 2006, 18 million metric tons were produced mostly by direct oxida- 
tion of ethylene, over 75 % of which were used for ethylene glycols production. 
Despite its technical and economical importance, experimental data on thcr- 
mophysical properties of ethylene oxide are rare, apart from basic properties 
at standard conditions [6]. This lack of data is mainly due to the hazardous 
nature of ethylene oxide. It is highly flammable, reactive, explosive at ele- 
vated temperatures, toxic, carcinogenic, and mutagenic. Therefore, it is an 
excellent example to show that molecular modeling and simulation can serve 
as a reliable route for obtaining thermophysical data in cases, where avoiding 
experiments is highly desirable. 

In the present work, a new molecular model for ethylene oxide was devel- 
oped. This model is based on prior work at our institute [7] and was further 
optimized to experimental vapor-liquid equilibria (VLB), i.e. saturated liq- 
uid density, vapor pressure, and enthalpy of vaporization. Using this model, 
phase equilibria, thermal, caloric, transport properties, and surface tension 
were predicted and compared to experimental results, where possible. 

2.1 Optimization of the molecular model for ethylene oxide 

The molecular model for ethylene oxide consists of three Lennard-Jones (LJ) 
sites (one for each methylene group and one for the oxygen atom) plus one 
static point dipole. Thus, the united-atom approach was used. Due to the fact 
that ethylene oxide is a small molecule, the internal degrees of freedom may 
be neglected and the model was assumed to be rigid. Experimental data on 
the molecular structure was used to specify the geometric location of the LJ 
sites, cf. [7]. The dipole was located in the center of mass. 

A set of five adjustable parameters, i.e., the four LJ parameters and the 
dipole moment, was optimized. This was done by a Newton scheme using 
correlations to experimental bubble density, vapor pressure, and enthalpy of 
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vaporization data over the full range of the VLE between triple point and 
critical point. 

Following the optimization procedure of StoU [7] , an optimized parameter 
set was obtained. The model parameters are listed in Table 1. In Figure 1, 
saturated densities, vapor pressure, and enthalpy of vaporization are shown. 
The present model describes the vapor pressure p„ , the saturated liquid den- 
sity p' , and the enthalpy of vaporization Ah^ with mean relative deviations 
in the complete VLE range from triple point to critical point of Spa = 1.5 %, 
Sp' = 0.4 %, and SAh^ = 1-8 %, respectively. 

Table 1. Coordinates and parameters of the present molecular model for ethylene 
oxide. The three Lennard- Jones sites are denoted by the molecular group which they 
represent, while the single dipolar site is denoted "Dipole". All coordinates are in 
principal axes with respect to the center of mass. The orientation of the electrostatic 
site is definded in standard Euler angles, where ip is the azimuthal angle with respect 
to the x-y-plane and 9 is the inclination angle with respect to the 2;-axis. 



Interaction 


X 


y 


2 a e/ks 


6 






Site 


A 


A 


A A K 


deg 


deg 


D 


CH2(1) 


0.78000 





-0.48431 3.5266 84.739 








CH2(2) 


-0.78000 





-0.48431 3.5266 84.739 








0(3) 








0.73569 3.0929 62.126 








Dipole 








— — 








2.459 




00 5.0 10.0 15.0 20.0 2OO 250 300 350 400 450 

p/moir'' 7/K 

Fig. 1. Vapor- liquid equilibrium of ethylene oxide. Saturated densities (left) and 
vapor pressure and enthalpy of vaporization (right): • simulation results, — exper- 
imental data [6]. 
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The model, as specified in Table 1, was used for calculating other properties 
discussed in the subsequent section. No further parameter adjustments were 
made and thus all results for properties except p' , and Ah^ are fully 
predictive. 

2.2 Results and comparison to experimental data 

Individual VLE simulation results from the present model between 235 and 
445 K deviate by less than 3 % in vapor pressure and heat of vaporization, 
and less than 0.5 % in saturated liquid density. 

Saturated densities, vapor pressure, and enthalpy of vaporization of ethy- 
lene oxide at 375 K are all reproduced by the present model with deviations 
below 1 % and therefore within the experimental error given in [6]. 

Critical properties are in excellent agreement with experimental data, be- 
ing within 1.3 %. The normal boiling temperature deviates by less than 0.5 % 
from experimental data for the present model. 

The second virial coefficient was calculated by direct evaluation of the 
intermolecular potential. It is underpredicted from simulation by 3.4 %. 

Experimental data on isobaric heat capacity and isothermal compressibil- 
ity for the saturated fluid states at 375 K have not been available to us prior 
to the deadline of the Simulation Challenge. Later, the IFPSC published ref- 
erence values, which are used for comparison here. In the saturated liquid, the 
simulation overpredicts isobaric heat capacity by 13 %, while in the saturated 
vapor, it is overpredicted by 18 %. Reference values for isothermal compress- 
ibility are obtained from the Brelvi-O'Connell correlation [8] for the saturated 
liquid and from the virial equation for the saturated vapor. Simulations for 
the present model yield a 39 % higher result for the liquid state and a 5.6 % 
higher result for the vapor state. 

While the deviation from the reference value for the isothermal compress- 
ibility in the saturated vapor state is within the assumed uncertainty of 9.2 %, 
the one in saturated liquid state is significantly above the assumed uncertainty 
of 23.1 %. This is astonishing as the densities predicted by the present model 
along the bubble line are in good agreement with experimental ethylene oxide 
data. An additional check on our value after the deadline of the Simulation 
Challenge, evaluating simulation runs at different pressures for 375 K, yielded 
a value within the statistical uncertainty of the previous simulation result. 
On the other side, the Brelvi-O'Connell correlation is known to give system- 
atic deviations to too low compressibilities for polar substances similar to 
ethylene oxide, e.g. sulfur dioxide or acetone, when using the critical density 
for reduction [8]. Therefore, the authors believe that molecular modeling and 
simulation yields more reliable results here than the standard method does. 

Experimental surface tension data are available between 200 and 296 K 
and may be extrapolated to higher temperatures using a correlation taken 
from [6]. The simulation results give deviations of —17 %. 
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Experimental transport properties for ethylene oxide are scarce. In fact, 
for viscosity and thermal conductivity at 375 K they are only available in 
the saturated vapor state. Here, the results from the present model agree 
with the experimental values within their statistical uncertainty. After the 
deadline of the Simulation Challenge, the IFPSC has given reference values for 
viscosity in the saturated liquid state from an extrapolation of experimental 
data from 223 to 282 K and for thermal conductivity from the Missenard 
method [9]. Simulation results agree with the reference values within their 
statistical uncertainties and the uncertainties of the reference values proposed 
by IFPSC. 

It can be seen that the molecular model for ethylene oxide is capable to rea- 
sonably predict the wide variety of properties. This underlines that molecular 
modeling and simulation in combination with the chosen molecular modeling 
route, using experimental bubble density, vapor pressure, and enthalpy of va- 
porization only, can be followed for predicting properties where experimental 
data is insufficient. 



3 An optimized molecular model for ammonia 

Ammonia is a well-known chemical intermediate, mostly used in fertilizer 
industries; another important application is its use as a refrigerant. Due to 
its simple symmetric structure and its strong intermolecular interactions it is 
also of high academic interest both experimentally and theoretically. 

In the present work a new molecular model for ammonia is proposed, which 
is based on the work of Kristof et al. [10] and improved by including data on 
geometry and electrostatics from ab initio quantum mechanical calculations. 

3.1 Selection of model type and parameterization 

For the present molecular model for ammonia, a single Lennard- Jones po- 
tential was assumed to describe the dispersive and repulsive interactions. The 
electrostatic interactions as well as hydrogen bonding were modeled by a total 
of four partial charges. This modeling approach was also followed by Kristof 
et al. [10] for ammonia. 

The Lennard-Jones site and the partial charges were placed according 
to the nucleus positions obtained from a quantum mechanical geometry op- 
timization. For an initial model, the magnitude of the charges were taken 
directly from quantum mechanics. This initial model was subsequently op- 
timized to experimental VLE data. The optimized parameters are given in 
Table 2. 

To optimize the molecular model, the two Lennard-Jones parameters were 
adjusted to experimental saturated liquid density, vapor pressure, and en- 
thalpy of vaporization using a Newton scheme as proposed by Stoll [7] . These 
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Table 2. Parameters of the molecular model for ammonia (e — 1.6021 -10 C is 
the elementary charge). 
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properties were chosen for the adjustment as they all represent major char- 
acteristics of the fluid region. Furthermore, they are relatively easy to be 
measured and are available for many components of technical interest. 

3.2 Results and discussion 
Vapor-liquid equilibria 

VLB results for the new model are compared to data obtained from a reference 
quality equation of state (EOS) [11] in Figure 2. This figure also includes the 
results that we calculated using the model from Kristof et al. [10]. 

The model of Kristof et al. shows noticeable deviations from experimental 
data. The mean unsigned errors over the range of VLE are 1.9 % in saturated 
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Fig. 2. Vapor-liquid equilibrium of ammonia. Saturated densities (left) and vapor 
pressure and enthalpy of vaporization (left): • simulation results from the new 
molecular model, ■ simulation results from the model by Kristof et al. [10], — 
experimental data [6]. 
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liquid density. 13 % in vapor pressure and 5.1 % in enthalpy of vaporization. 
With the new model, a significant improvement was achieved compared to 
the model from Kristof et al. The description of the experimental VLE is very 
good, the mean unsigned deviations in saturated liquid density, vapor pressure 
and enthalpy of vaporization are 0.7, 1.6, and 2.7 %, respectively. 

Mathews [12] gives experimental critical values of temperature, density and 
pressure for ammonia: Tc =405.65 K, pc =13.8 mol/1, and pc =11.28 MPa. 
Following the procedure suggested by Lotfi et al. [13] the critical properties 
Tc =395.82 K, pc =14.0 mol/1, and Pc =11.26 MPa for the model of Kristof et 
al. were calculated, where the critical temperature is underestimated by 2.4 %. 
For the new model Tc =402.21 K, =13.4 mol/1, and Pc =10.52 MPa were 
obtained. The new model gives very good results for the critical temperature, 
while it underpredicts the critical pressure slightly. 

Homogeneous region 

In many technical applications thermodynamic properties in the homogeneous 
fluid region are needed. Thus, the new molecular model for ammonia was 
tested on its predictive capabilities for such states. 

Thermal and caloric properties were predicted with the new model in the 
homogenous liquid, vapor and supercritical fluid region. In total, 70 state 
points were regarded, covering a large range of states with temperatures of 
up to 700 K and pressures of up to 700 MPa. In Figure 3, relative deviations 
between simulation and reference EOS [11] in terms of density and enthalpy 
are shown. The deviations are typically below 3 % for density with the excep- 
tion of the extended critical region, where a maximum deviation of 6.8 % is 
found, and below 5 % for enthalpy. 

These results confirm the modeling procedure. By adjustment to VLE data 
only, quantitatively correct predictions in most of the technically important 
fluid region can be obtained. 

Structural quantities 

Due to its scientific and technical importance, experimental data on the mi- 
croscopic structure of liquid ammonia are available. Ricci et al. [14] applied 
neutron diffraction and published all three types of atom-atom pair correla- 
tions, namely nitrogen- nitrogen (N-N), nitrogen-hydrogen (N-H), and hydro- 
gen-hydrogen (H-H). In Figure 4, these experimental radial distribution func- 
tions for liquid ammonia at 273.15 K and 0.483 MPa are compared to present 
predictive simulation data based on the new ammonia model. 

It was found that the structural properties are in very good agreement, 
although no adjustment was done regarding these properties. The atom-atom 
distance of the first three layers is predicted correctly, while only minor over- 
shootings in the first peak are found. Please note that the first experimental 
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Fig. 3. Relative deviations for the density (left) and the enthalpy (right) between 
simulation and reference EOS [11] ijiz = (zsim — 2eos)/2cos) in the homogeneous 
region: O simulation data of new model, — vapor pressure curve. The size of the 
bubbles denotes the relative deviation as indicated in the plot. 




Fig. 4. Pair correlation functions of ammonia: — simulation data with the new 
model, o experimental data [14]. 
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peak in (?n-h ^^nd gn-H show intramolecular pair correlations, which are not 
included in the present model. 

In the experimental radial distribution function .gN-H the hydrogen bond- 
ing of ammonia can be seen at 2-2.5 A. Due to the simplified approximation 
by eccentric partial charges, the molecular model is not capable to describe 
this effect completely. But even with this simple model a small shoulder at 
2.5 A is obtained. 



4 Nucleation processes during condensation 

Understanding homogeneous nucleation is required to develop an accurate 
theoretical approach to nucleation that extends to more complex and tech- 
nically more relevant heterogeneous systems. MD simulations can be used to 
investigate the condensation of homogeneous vapors at high supersaturations. 

The nucleation rate J is influenced to a large extent by the surface energy 
of droplets, which also determines how many droplets are formed and from 
which size on they become stable. The accuracy of different theoretical ex- 
pressions for the surface energy can be assessed by comparison to results from 
MD simulations. 

4.1 Nucleation theory 

Classical nucleation theory (CNT) was developed by Volmer and Weber [15] in 
the 1920s and further extended through many contributions over the following 
decades [16]. It is founded on the capillarity approximation: droplets emerging 
during nucleation are assumed to have the same thermodynamic properties 
as the bulk saturated liquid. In particular, the specific surface energy of the 
emerging nano-scaled droplets is assumed to be the surface tension of the 
planar phase boundary in equilibrium. 

Laaksonen, Ford, and Kulmala [17] (LFK) proposed a surface energy coef- 
ficient that depends on the number of molecules in the droplet, i.e. the size of 
the droplet. Tanaka et al. [18] found that this expression leads to nucleation 
rates which agree with their simulation results. It was shown both theoret- 
ically [19, 20] and by simulation [21, 22] that the surface tension acting in 
the curved interface of nano-scaled droplets is actually much lower than in a 
planar interface. This finding underlines the necessity of a correction to CNT. 

4.2 Simulation method 

Both the critical droplet size and the nucleation rate can be determined by 
molecular simulation. After an initial period of equilibration, the steady state 
distribution of droplets is established. The critical droplet size, i.e. the size 
from which on droplets more probably grow than decay, can either be calcu- 
lated from the equilibrium distribution [23] or by the "nucleation theorem" 
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[24]. Both methods require data on the nucleation rate J, i.e. the number 
of macroscopic droplets emerging per volume and time in a steady state at 
constant supersaturation. From an MD simulation of a supersaturated vapor, 
this rate can straightforwardly be extracted by counting the droplets that 
exceed a certain threshold size. This method was proposed by Yasuoka and 
Matsumoto [23], who found that as long as this threshold is larger than the 
critical droplet size, its precise choice hardly affects the observed value of J . 

The nucleation rate J is the number of nuclei formed per volume and time 
in a supersaturated vapor 

^(') = E^' (1) 

for a threshold size t larger than the critical size, measured after an initial 
temporal delay required to form sufficiently large nuclei. For nucleation rates 
obtained from molecular simulations, a very low size threshold is used, roughly 
of the same order of magnitude as the critical size. Because the choice of the 
threshold may influence the result, it must be indicated explicitly, i.e. as J(t) 
instead of J. 

For simulations in the canonical ensemble, it has to be taken into ac- 
count that as the condensation proceeds, the density of the remaining vapor 
decreases and the pressure in the vapor is reduced significantly. This causes 
larger nuclei to be formed at a lower rate. Thus, the nucleation rates are given 
together with pressure values which were taken in the middle of the interval 
where the value of J(t) was obtained by linear approximation. 

For the present study, methane was modeled as a LJ fluid and ethane as 
a rigid two-center L J fluid with an embedded point quadrupole [25] . Further- 




Fig. 5. Nucleation rates of methane and ethane at low temperatures from simulation 
according to the Yasuoka-Matsumoto method and different threshold sizes (■) as 
weU as CNT (— ) and LFK (---). 
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Fig. 6. Nucleation rate of the LJTS fluid over supersaturated pressure from tlie 
present direct simulations for different threshold values (• i = 25, O i = 50, ■ i G 

{75, 100}, A i > 150) and following the standard CNT ( ), the LFK modification 

of CNT (■ ■ ■), and the new SPC modification of CNT (— ). 

more, the truncated and shifted LJ (LJTS) fluid with a cutoff radius of 2.5 
times the size parameter a was studied. The LJTS fluid is an accurate model 
for fluid noble gases and methane [21]. 

Within a typical computation time of about 24 hours, a time interval of 
a few nanoseconds can be simulated for systems with a volume of 10^^^ m'^. 
Given that the nucleation rate is determined as units per volume and time, 
only values that exceed 10"^" m~'^s~^ can be obtained from molecular simula- 
tion. However, only nucleation rates below 10^^ m~'^s~^ can be measured in 
experiments at present. Therefore, computational power is crucial to keep the 
gap between simulation and experiment as small as possible. 

In Figure 5 the nucleation rates from MD simulation are shown in com- 
parison to predictions from CNT and LFK theories. At the given state points, 
simulation results confirm both CNT and LFK, deviations are throughout 
lower than three orders of magnitude. 

4.3 Simulation results 

However, from a series of simulations it was seen that standard CNT only 
predicts the nucleation rate J with an acceptable accuracy, but leads to de- 
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viations for the critical droplet size; the LFK modification provides excellent 
predictions for the critical size but not for the temperature dependence of J. 
LFK underpredicts the nucleation rate at high temperatures by several orders 
of magnitude. Both theories assume an inappropriate curvature dependence 
of the surface tension, although for this essential property of inhomogcncous 
systems a qualitatively correct expression is known since the 1940s [26]. With 
the collected simulation data on critical droplet size and nucleation rate over 
a broad range of temperatures, enough quantitative information is available 
to formulate a more adequate modification of CNT. For the LJTS fiuid, a 
surface property corrected (SPC) modification of CNT was developed which 
postulates a non-spherical surface and is based on surface tension values ob- 
tained from molecular simulation. As shown in Figure 6, the SPC modification 
accurately describes the nucleation rate of the LJTS fiuid. 

5 Computing performance 

All simulations presented in Sections 2.1 and 3 were carried out with the 
MPI based molecular simulation program ms2 developed in our group. The 
parallelization of the molecular dynamics part of ms2 is based on Plimptons 
force decomposition algorithm [27]. The nucletion simulations were carried 
out with the massivly parallel program Isl. Here, the parallelization relies on 
spatial domain decomposition for parallel simulations [28] due to the large 
size of the simulated systems (up to 2 million particles). 

With ms2 typical simulation runs to determine the vapor-liquid equilib- 
rium employ 4-8 CPUs running for 4-6 hours. For model optimization or the 
comprehensive study on mixtures, a large number of independent simulations 
are necessary and can be performed in parallel. For the prediction of other 
thermophysical properties, additional simulations were performed using up to 
32 CPUs and running up to 72 hours depending on the system size and desired 
property. 

For the simulation of nucleation processes, very large systems must be con- 
sidered. Here, the massively parallel program Isl was used. Table 3 demon- 
Table 3. Scaling of the massively parallel program Isl on HP XC4000. Particle 
numbers A*' are 0.5-2 million. 
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strates the good scaling of the program on the HP XC4000 cluster at the 
Steinbuch Center for Computing, Karlsruhe. 
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